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The properties of water can have a strong dependence on the confinement. Here, we consider a water 
monolayer nanoconfined between hydrophobic parallel walls under conditions that prevent its 
crystallization. We investigate, by simulations of a many-body coarse-grained water model, how the 
properties of the liquid are affected by the confinement. We show, by studying the response functions and 
the correlation length and by performing finite-size scaling of the appropriate order parameter, that at low 
temperature the monolayer undergoes a liquid-liquid phase transition ending in a critical point in the 
universality class of the two-dimensional (2D) Ising model. Surprisingly, by reducing the linear size L of the 
walls, keeping the walls separation h constant, we find a 2D-3D crossover for the universality class of the 
liquid-liquid critical point for L/h ~ 50, i.e. for a monolayer thickness that is small compared to its extension. 
This result is drastically different from what is reported for simple liquids, where the crossover occurs for 
L//i~5, and is consistent with experimental results and atomistic simulations. We shed light on these 
findings showing that they are a consequence of the strong cooperativity and the low coordination number 
of the hydrogen bond network that characterizes water. 

The study of nanoconfined water is of great interest for applications in nanotechnology and nanoscience 1 . The 
confinement of water in quasi-one or two dimensions (2D) is leading to the discovery of new and contro- 
versial phenomena in experiments 15 and simulations 4,6,7 . Nanoconfinement, both in hydrophilic and hydro- 
phobic materials, can keep water in the liquid phase at temperatures as low as 130 K at ambient pressure 2 . At these 
temperatures T and pressures P experiments cannot probe liquid water in the bulk, because water freezes faster 
then the minimum observation time of usual techniques, resulting in an experimental "no man's land" 8 . 
Nevertheless, new kind of experiments 910 and numerical simulations 1 1 can access this region, revealing interesting 
phenomena in the metastable state. In particular, Poole et al. found, by molecular dynamics simulations of 
supercooled water, a liquid-liquid critical point (LLCP), in the "no mans land", at the end of a first-order 
liquid-liquid phase transition (LLPT) line between two metastable liquids phases with different density p: the 
high-density liquid (HDL) at higher TandP, and the low-density liquid (LDL) at lower Tand P u . The presence of 
a LLPT is experimentally observed in other systems 12-21 , consistent with theoretical models fitted to water 
experimental data 22-24 , and is recovered by simulations of a number of models of water 11,25 31 and other anomalous 
liquids 32-37 . Alternative ideas, and their differences, have been discussed 38-42 , and it has been debated if experi- 
ments on confined water in the "no man's land" can be a way to test these ideas 2 , motivating several theoretical 
works 43 . 

Here, to analyze the thermodynamic properties of water in confinement we consider a water monolayer 
between hydrophobic walls of area L 2 separated by h ~ 0.5 nm (Fig. 1). Atomistic simulations 7 show that water 
under these conditions does not crystallize, but arranges in a disordered liquid layer, whose projection on one of 
the surfaces has square symmetry, with each water molecule having four nearest neighbors (n.n.). The molecules 
maximize their intermolecular distance by adjusting at different heights with respect to the two walls. 

We adopt a many-body model that reproduces water properties 31,40,44-50 . We simulate ~ 10 5 state points, each 
with statistics of 5 X 10 6 independent calculations, for systems with N = 2.5 X 10 3 , . . ., 1.6 X 10 5 water molecules 
at constant N, P and T, using a cluster Monte Carlo algorithm 46-48 , for a wide range of T and P. All quantities are 
calculated in internal units, as described in the Methods section. 

Results 

We calculate the density p = N/Vof the system as function of T along isobars. For a broad range of P, we find a 
maximum and a minimum of density along each isobar (Fig. 2a) according to experimental evidences for bulk and 
confined water 52 . These maxima and minima identify, for each P, the temperature of maximum density (TMD) 
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Figure 1 | Schematic view of a section of the water monolayer confined 
between hydrophobic walls of size L x L separated by h ~ 0.5 nm. 

and the temperature of minimum density (TminD). The TMD locus 
merges the TminD line as in experiments 52 and other models 53 . 

At low T a discontinuous change in p is observed for 
1 >Pvo/(4e) >0.5, where the parameters v 0 and £ are explained in 
the Methods section, as it would be expected in correspondence of 
the HDL-LDL phase transition. At very high pressures 
(Pv 0 /(4e) > 1) the system behaves as a normal liquid, with mono- 
tonically increase of p upon decrease of T. 

We estimate the liquid-to-gas (LG) spinodal at Pv 0 /(4e) <0 for 
low T (Fig. 2) as the temperature along an isobar at which we find a 
discontinuous jump of p to zero value by heating the system. The LG 
spinodal identifies the locus of the stability limit of liquid phase with 
respect to the gas phase: at pressures below the LG spinodal in the 
P — T plane is no longer possible to equilibrate the system in the 
liquid phase. The LG spinodal continues at positive pressures ending 
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Figure 2 | (a) Isobaric density variation for 10" water molecule. Lines join 
simulated state points (—150 for each isobar). P increases from —0.5 
(bottom curve) to 1.5 (4e) /vn (top curve). Along each isobar we locate the 
maximum p (green squares at high T) and the minimum p (green small 
circles at low T) and the liquid-gas spinodal (open large circles at low P). 
(b) Loci of TMD, TminD, liquid-gas spinodal and liquid-liquid spinodal in 
T — P plane. Dashed lines with labels represent the isochores of the system 
from pv 0 = 0.43 (bottom) to pv 0 = 0.80 (top). Dashed lines without labels 
represent intermediate isochores. TMD and TminD correspond to the loci 
of minima and maxima, respectively, along isochores in the T — P plane. 
We estimate the critical isochore at pv 0 ~ 0.47 (red circles). All the 
isochores with 0.47 < pv 0 < 0.76 intersect with the critical isochore for 
Pvo/(4e) >0.5 along the LL spinodal (tick turquoise) line. 



in the LG critical point (data not shown). We observe that the TMD 
line approaches the LG spinodal, without touching it (Fig. 2). We 
recover the LG spinodal also as envelope of isochores (Fig. 2b). 

We find a second envelope of isochores at lower T and higher P, 
pointing out the liquid-to-liquid (LL) spinodal. Indeed, the two spi- 
nodals associated to the LLPT, i.e. the HDL-to-LDL spinodal and the 
LDL-to-HDL spinodal, collapse one on top of the other and are 
indistinguishable within our numerical resolution. Nevertheless, 
we clearly see that isochores are gathering around the points 
(Tfc 5 /(4£) -0.06, Pv 0 /(4e) -0.5) and (Tfc B /(4e)=0, Pv 0 /(4£) = 1), 
where k B is the Boltzmann constant, marking two possible critical 
regions (Fig. 2b). 

We calculate the isothermal compressibility by its definition K T = 
— (1/(V)) (d(V)/8P) T and by the fluctuation-dissipation theorem K T 
= (AV 2 )/k B TV along isobars, Kj(T), and along isotherms, Kj(P) 
(Fig. 3), where (V) = V is the average volume and (AV 2 ) the volume 
fluctuations. We find two loci of extrema for each quantity Kj( T) and 
K T (P): one of strong maxima and one of weak maxima. The loci of 
strong maxima in Kj(T) and K r (P), respectively rv^ Max (T) and 
jq, Max (P), overlap within the error bar with the LL spinodal. The 
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v T yi j auva iv r y± j increase m me range 
Pv 0 /(4e)e[0.55;0.6] and Tfc B /(4£)e[0.05; 0.06] (Fig. 3), consistent 
with the existence of a critical region. The stronger maxima disappear 
forPv 0 /(4£)<0.4. 

We find also loci of weak maxima, K™ Max (T) and K™ Max (P) and 
minima K™ ln ( T) and K™ m (P) . The loci of weak extrema and minima 
of Kj(T) and Kj{P) do not coincide in the T — P plane. The locus of 
weak maxima along isotherms K* Max (P) merges with the locus of 
minima rC™ m (P) at the point where the slope of both loci is dPldT^> 
oo. Furthermore, both loci approach to the LL spinodal at high P. The 
locus of weak maxima along isobars K™ Max (T) approaches the LL 
spinodal where K T exhibits the strongest maxima, and merges with 
the locus of minima rC™ ln (T) where the slope of both loci is dPlcT^> 
0 (data at high P and Tnot shown in Fig 3). This locus intersects the 
TMD at its turning point. Indeed, as reported in Ref. 39 and in the 
Methods section, the temperature derivative of isobaric K T along 
the TMD line is related to the slope of TMD line 



8Kj 



1 (d 2 V/dT 2 ] 
V {dP/dT\ 



(1) 



where all the quantities are calculated along the TMD line. Hence the 
locus of extrema in Kj{T), where (dK T /dT) P = 0, crosses the TMD 
line where the slope (SP/SPJtmd is infinite. We observe also that the 
weak maxima of Kj{T) and K T (P) increase as they approach the LL 
spinodal. All loci of extrema in K T are summarized in Fig. 3. 

Next we calculate the isobaric specific heat C P = (6(H)/ dT) P = 
(AH 2 )/k B T along isotherms and isobars, where (H) = (j4f) +P{V) is 
the average enthalpy, .yf is the Hamiltonian as defined in the 
Methods section, (AH 2 ) is the enthalpy fluctuations (Fig. 4). We find 
two maxima at low P separated by a minimum. At high- T the maxima 
are broader and weaker than those at low-T. As discussed in Ref. 49, 
the maxima at high T are related to maxima in fluctuations of the HB 
number N HB , while the maxima at low T are a consequence of max- 
ima in fluctuations of the number N coop of cooperative HBs. The lines 
of strong maxima at constant P and constant T, respectively Cp Max (T) 
and Cp Max (P), overlap for all the considered pressures, and both 
maxima are more pronounced in the range Pv 0 /(4£)e[0.5, 0.6] and 
7V(4£)e[0.06, 0.07]. The weak maxima C™ Max (P) and C£ Max (T) 
increase approaching the LL spinodal and have their larger maxima at 
the state point where they converge to the strong maxima, consistent 
with the occurrence of a critical point for a finite system (Fig. 4). The 
lines of weak maxima overlap for all positive pressures, branching off 
at negative pressures. At negative pressures, the locus Cp Max (P) bends 
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Figure 3 | (a) Loci of strong maxima (i"Cy Max (T)), weak maxima (if™ Max (T) in the inset) and minima (i<C™ ln (T) marked with large triangles in the inset) 
along isobars for Kj(T). (b) Loci of strong maxima (K^ Max (P)), weak maxima (if™ Max (P) in the inset) and minima (iC™ m (P) marked with large triangles in 
the inset) along isotherms. The weak maxima merge with minima, (c) Projection of extrema of K T m T — P plane. The strong maxima (symbols), 
weak maxima (solid lines) and minima (dashed lines) of Kj( T) (orange) and Kj{P) (blue) form loci in T — P plane that relate to each other and intersect 
with the TMD line following the thermodynamic relations discussed in the text. The large yellow circle with label A identifies the region where i<Cy Max (T) 
and K"y Max (P) converge and display the largest maxima, consistent with the occurrence of a critical point in a finite-size system. Symbols not explained 
here are as in Fig. 2. 



toward the turning point of the TMD line, as discussed in Methods 
section and in Ref. 53. Indeed, according to the relation 



dC P 

~8P 



-t(£\ (£L) 



\dTj TMD \8PdTj TMD 



(2) 



in case of intersection between the locus of extrema (dCp/dP) T = 0 
and the TMD line, it results that (dP/8T) TMD = 0. Note that, as we 
explain in the Methods section, the relation (2) does not imply any 
change in the slope of the TminD line at the intersection with the 
locus of {dC P /dP) T = 0. 

We calculate also the thermal expansivity a P = (l/( V)) (d(V)/dT) P 
along isotherms and isobars (Fig. 5). As for the other response func- 
tions, we find two loci of strong extrema, minima in this case, 
« p mm (P) and ap mm (T), along isotherms and isobars, respectively 
showing a divergent behavior in the same region where we find the 
strong maxima of K T and Cp. From this region two loci of weaker 
minima depart. We find that the locus of weak minima along isobars 
a wmin (• r -j b enc j s towarc j t he turning point of the TMD. Although our 
calculations for a P do not allow us to observe the crossing with the 
TMD line, based on the relation (see Methods) 



(—) 
tmd \dPSTj TMD 



(3) 



that holds at the TMD line, we can conclude that ap mn (r) should 
have zero T-derivative if it crosses the point where the TMD turns 



into the TminD line, because in this point the TMD slope approaches 
zero. 

The locus of weaker minima along isotherms <y™ m (P), merges 
with the locus of maxima a^ ax (P) at the state point where the slope 
of both loci is dP/dT— > 00 (not shown in Fig. 5). According to the 
thermodynamic relation, discussed in Methods section, 

8oip\ f dK T ^ 

~8P L = ~ [~dT 



(4) 



we find that the locus of extrema in thermal expansivity along iso- 
therms coincides, within the error bars, with the locus of extrema of 
isothermal compressibility along isobars (Fig. 5c). 

All the loci of extrema of response functions that converge toward 
the same region A in Fig. 3, 4 and 5 increase in their absolute values. 
Because the increase of response functions is related to the increase of 
fluctuations and this is, in turn, related to the increase of correlation 
length f, to estimate C we calculate the spatial correlation function 

[( CT #;) ff flt(n))-( ff 9 



G(r) 



4N„ 



'i-nl= r 



5) 



where r, is the position of the molecule i, \% — ?i\=T the distance 
between molecule i and molecule / and (•) the thermodynamic aver- 
age. The states of the water molecule, as well as the density p, the 
energy E and the entropy S of the system, are completely described 
by the bonding variables a^. Therefore, the function G(r) accounts 
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Figure 4 | (a) Loci of strong maxima (Cp Max (T)) and weak maxima (Cp Max (T) in the inset) along isobars for C P . (b) Loci of strong maxima (Cp Max (P)) and 
weak maxima (C™ Max (P) in the inset) along isotherms, (c) Projection of C P maxima in T — P plane. The large circle with A identifies the region 
where C P shows the strongest maximum. Symbols not explained here are as in Fig. 2. 



for the fluctuations in p, E and S and allows us to evaluate the 
correlation length because the order parameter of the LLPT, as we 
discuss in the following, is related to a linear combination of p and E. 
Note that, instead, the density-density correlation function would 
give only an approximate estimate of {. 

We observe an exponential decay of G(r) ~ e~ rli at high tempera- 
tures in a broad range of pressures. Approaching the region A, the 
correlation function can be written as G(r) ~ e _rf{ /r rf ~ 2+ '' where d is 
the dimension of the system and r\ a (critical) positive exponent. 
When £, is of the order of the system size, the exponential factor 
approaches a constant leaving the power-law as the dominant con- 
tribution for the decay. 

At P below the region A, we find that t, has a maximum, £ Max , along 
isobars and that i; Max increases approaching A (Fig. 6). The C Ma * locus 
coincides with the locus of strong extrema of Cp, K T and a P (Fig. 6b). 
We observe that this common locus converges to A and that all the 
extrema increase approaching A. This behavior is consistent with the 
identification of A with the critical region of the LLCP. Furthermore, 
we identify the common locus with the Widom line that, by defini- 
tion, is the C Max locus departing from the LLCP in the one-phase 
region 54,55 . Our calculations allow us to locate the Widom line at 
any P down to the liquid-to-gas spinodal. 

At P above the region A, we find the continuation of the c Max line, 
but with maxima that decrease for increasing P, as expected at the LL 
spinodal that ends in the LLCP (Fig. 6). Therefore, we identify the 
high-P part of the C Max locus with the LL spinodal. Along this line the 
density, the energy and the entropy of the liquid are discontinuous, as 
discussed in previous works 31,40,44-49 . 

To better locate and characterize the LLCP in A we need to define 
the correct order parameter (o.p.) describing the LLPT. According to 



mixed-field finite-size scaling theory 56 , a density-driven fluid-fluid 
phase transition is described by an o.p. M = p* + su*, where p* 
= pv 0 represents the leading term (number density), u = E/(eN) is 
the energy density (both quantities are dimensionless) and s is the 
mixed-field parameter. Such linear combination is necessary in order 
to get the right symmetry of the o.p. distribution Gjv(M) at the critical 
point where Q N (M) ocp d (x). Here is x = B(M - M c ), B = a^N ll/d \ 
P is the critical exponent that governs M, v is the critical exponent that 
governs with v and /i defined by the universality class, a M is a non- 
universal system-dependent parameter and pd is an universal function 
characteristic of the Ising fixed-point in d dimensions. We adjust B 
and M c so that Qn(M) has zero mean and unit variance. 

We combine, using the multiple histogram reweighting method 57 
described in the Methods section, a set of 3 X 10 4 MC independent 
configurations for ~ 300 state points with 0.040< Tk B /{^t) <0.065 
and 0.40 < Pv 0 /(4£) < 0.75. We verify, by tuning s, T and P, that there 
is a point within the region A where the calculated QjvOc) has a sym- 
metric shape with respect to x = 0 (Fig. 7). We find s = 0.25 ± 0.03 for 
our range of N. The resulting critical parameters T C (N), P C (N) and the 
normalization factor B(N) follow the expected finite-size behaviors 
with 2D Ising critical exponents 56 . From the finite-size analysis we 
extract the asymptotic values T c k B /(4e) =0.0597 + 0.0001 and 
P c v 0 /(4£) =0.555 + 0.002. 

The presence of a first order phase transition ending in a critical 
point, associated to the o.p. M, is confirmed by the finite size analysis 
of the Challa-Landau-Binder parameter 58 of M 



U A 



(6) 
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where the symbol refers to the thermodynamic average for a 
system with N water molecules. U M quantifies the bimodality in 
Q;v(M). The isobaric value of U M shows a minimum at the temper- 
ature where Q;v(M) mostly deviates with respect to a symmetric 
distribution (Fig. 8). Minimum of U M converges to 2/3 in the ther- 
modynamic limit away from a first order phase transition, while it 
approaches to a value <2/3 where the bimodality of Q N (M) indicates 
the presence of phase coexistence. 

These results are consistent with the behavior of the Gibbs free 
energy G calculated with the histogram reweighting method (Fig. 9). 
In particular, we calculate G along isotherms, for P crossing the LLPT 
and the loci of weak maxima in Kj{T) and C P (P). We find that the 
behavior of G for T < T c is consistent with the occurrence of a 
discontinuity in volume V = dG/dP, in the thermodynamic limit, 
with a decrease of V corresponding to the transition from LDL to 
HDL for increasing P. Crossing the loci Kr(r) wMax and Cp(P) wMax 
the volume decreases with pressure without any discontinuity as 
expected in the one-phase region. 

The distribution Qn(N) adjust well to the data only for large N. 
We, therefore, perform a more systematic analysis. For each N, we 
quantify the deviation of the calculatedp(N) from the expected pi for 
the 2D Ising. Furthermore, due to the behavior of data for small N 
(Fig. 7a), we calculate the deviation from the 3D Ising pi 56 . We 
estimate the Kullback-Leibler divergence 51,59 , 



of the probability distribution pi{N) of x, from the theoretical value 
pdj of Xj (i = 1, ...,«) in d dimensions (Fig. 10a), and the Liu et al. 
deviation 51 , 

WdWs iu«2/mmzM (8 ) 

n pd, v eak—pd,x=0 

withp^peak —pd,x=o difference between the distribution peak and its 
value at x = 0 (Fig. 10b). 

We confirm 5 ~ 0.25 forp 2 and find s = 0.10 ± 0.02 for p 3 for our 
range of N. For both D|p" and Wj, with d = 2 and d = 3, we find 
minima at T c ks/(4e) ~0.06 andP c vo/(4£) ~0.55 that become stron- 
ger for increasing N. We find that and W 2 decrease with increas- 
ing N, vanishing for N — » 00 (Fig. 10). Therefore, for an infinite 
monolayer between hydrophobic walls separated by h ~ 0.5 nm, 
the system has a LLCP that belongs to the 2D Ising universality class, 
as expected from our representation of the system as the 2D projec- 
tion of the monolayer. 

However, by increasing the confinement, i.e. reducing N and L at 
constant p, Df L and W 2 become larger than and W 3 , respect- 
ively. Therefore, the calculatedp(N) deviates from the 3D probability 
distribution less than from the 2D probability distribution. For N = 
2500 we find that both Df h and W 3 have values approximately equal 
to those for Df h and W 2 calculated for a system ten times larger. In 
particular we find Df-~0 for N = 2500. Hence, by increasing the 
confinement of the monolayer at constant p, the LLCP has a behavior 
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Figure 6 | (a) The correlation length t along isobars for N = 10" water 
molecules has maxima that increase for P approaching the critical region A. 
(b) The locus of £ maxima coincides with the loci of strong extrema of K T , 
C P and n P . The Widom line is by definition the locus of c maxima at high T 
departing from the LLCP, that we locate within the critical region A, as 
discussed in the text. 

that approximates better the bulk 25-3038 , with a crossover between 2D 
and 3D-behavior occurring at IV ~ 10 4 . 

This dimensional crossover is confirmed by the finite-size analysis 
of the Gibbs free energy cost AG/(k B T c ) to form an interface between 
the two liquids in the vicinity of the LLCP, calculated as 
AG(N)=-k B T c (N)[ln^(j^V)-]n^(j^,V)], where gff 
and are the minimum and maximum values of the probability 

distribution .^^{j^ \V) of configurations of N water molecules with 
energy Jif and volume V at the LLCP. This quantity is expected to 
scale as AGccAP 1 ". We find that our data can be fitted as N* for small 
sizes and as for large sizes with a crossover around N = 10 4 
(Fig. 10c). Considering the value of the estimated p c in real units 
(~ lg/cm 3 ) 45 , the corresponding crossover wall-size is L~25 nm. 




-2-1012 

Rescaled order parameter x=B(M-M ) 




ln(A0 



Figure 7 | (a) The size-dependent probability distribution Q N for the 
rescaled o.p. x, calculated for T C (N), P C (N) and B(N), has a symmetric 
shape that approaches continuously (from N = 2500, symbols at the top at 
x = 0, to N = 40000, symbols at the bottom) the limiting form for the 2D 
Ising universality class (full blue line) and differs from the 3D Ising 
universality class case (full black line). Error bars are smaller than the 
symbols size, (b) The size-dependent LLCP temperature T C (N) and (c) 
pressure P C (N) (symbols), resulting from our best-fit of Q N , extrapolate to 
T c k B /(4€) ~ 0.0597 and P c v 0 /(4e) ~ 0.555, respectively, following the 
expected linear behaviors (lines), (d) The normalization factor B{N) 
(symbols) follows the power law function (dashed line) oc W"*. We use 
the d = 2 Ising critical exponents: 8 = 2 (correction to scaling), v = 1 and 
P = 1/8 (both defined in the text). 



Discussion 

Our rationale for this dimensional crossover at fixed h is that, when 
LI h decreases toward 1 , the characteristic way the critical fluctuations 
spread over the system, i.e. the universality class of the LLCP, resem- 
bles closely the bulk because the asymmetry among the three spatial 
dimensions is reduced. A similar result was found recently by Liu et 
al. for the gas-liquid critical point of a Lennard-Jones (LJ) system 
confined between walls by fixing L and varying h 51 . However, in the 
case considered by Liu et al. the crossover was expected because the 
number of layers of particles was increased from one to several, 
making the system more similar to the isotropic 3D case. Here, 
instead, we consider always one single layer, changing the proportion 
L/h by varying L. Therefore, it could be expected that the system 
belongs to the 2D universality class for any L. 

Furthermore, the extrapolation of the results for the LJ liquid to 
our case of a monolayer with h/ro — 1 .7, where r 0 is the water van der 



Waals diameter, would predict a dimensional crossover at L/h m 5 51 . 
Here, instead, we find the crossover at L/h st 50, i.e. one order of 
magnitude larger than the LJ case. We ascribe this enhancement of 
the crossover to (i) the presence of a cooperative HB network and (ii) 
the low coordination number that water has in both the monolayer 
and the bulk. These are the main differences between water and a LJ 
fluid. The cooperativity intensifies drastically the spreading of the 
critical fluctuations along a network, contributing to the effective 
dimensionality increase of the confined monolayer. Moreover, the 
HB network has in 3D a coordination number (z = 4) as low as in 2D, 
making the first coordination shell similar in both dimensions. 

Our findings are consistent with recent atomistic simulations of 
water nanoconfined between surfaces. 60 62 . Zhang et al. found that 
water dipolar fluctuations are enhanced in the direction parallel to 
the confining surfaces (hydrophobic graphene sheets) within a dis- 
tance of 0.5 nm 60 . Ballenegger and Hansen found similar results for 
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Figure 8 | Challa-Landau-Binder parameter U M (defined in the text) of the o.p. M for different system sizes, calculated for three pressures: (a) 
Pv 0 /(4c) = 0.9, (b) Pv 0 1 (4c) = 0.7, and (c) Pv 0 j (4c) =0.5 slightly below P c v 0 / (4c) ~ 0.555. The curves are calculated with the histogram reweighting 
method, (d) Scaling of the minima of U M for different P. The arrow points to value 2/3 corresponding to the absence of a first-order phase transition in the 
thermodynamic limit. Error bars are calculated propagating the statistical error from histogram reweighting method. 
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Figure 9 | Gibbs free energy G along isotherms, as function of P. Points are shifted so that G = 0 at the lowest P. Lines are guides for the eyes, (a) For 
T = 0.02 (4c I kg) < T c there is a discontinuity in the P-derivative of Gat P~ 0.952 (4c / Vo) >P C as expected at the LLPT, consistent with the behavior of 
the response functions at this state point (e.g., in Fig. 3b, 4b). (b) For T = 0.04 (4c/kB)< T c we observe the discontinuity in the P-derivative at 
P~ 0.865 (4c /vo) >P C , again consistent with the LLPT. The LDL has a lower chemical potential (/i = GIN) than the HDL, j.i LDL < Hhdl> due to the HB 
energy gain in the LDL. For T = 0.061 (4c /kg) (c) and for T = 0.15 (4c /ks) (d), both larger than T a we instead do not observe any discontinuity in the P- 
derivative of G by crossing the locus of Cp(P) wMax and the locus of Kr(T) wMaK , respectively, as expected in the one-phase region. 



SCIENTIFIC REPORTS | 4 : 4440 | DOI: 1 0. 1 038/srep04440 



7 




-4 



l.OxlO 4 1,0x10 s 

Number of water molecules N 



Figure 10 | (a) Kullback-Leibler divergence D^ L (N) and (b) Liu et al. deviations of the calculated p(N) from the Ising universal function^ in d = 2 
(open symbols) and d = 3 (closed symbols), as a function of 1/iV, with N water molecules, at constant p^p c . In both panels lines are power-law fits 
and we observe a crossover between 2D and 3D behavior at N ~ 10 4 . (c) The free-energy cost to form an interface between the two liquids coexisting at the 
LLCP scales as AGocN^ with d = 3 for N < 10" and d = 2 for AT > 10". 



confined polar fluids, including water, within ~ 0.5 nm distance 
from the hydrophobic surface 61 . Bonthuis et al. extended these results 
to both hydrophilic and hydrophobic confining surfaces. All these 
findings are consistent with our result showing the enhancement of 
the fluctuations of the o.p. in the direction parallel to the confining 
walls separated by h ~ 0.5 nm. Furthermore, Zhang et al. observed 
that the effect does not depend on the details of the water-surface 
interaction but stems from the very presence of interfaces 60 . This is 
confirmed by our study, where the water-interface interaction is 
purely due to excluded volume. Following the authors of Ref. 60, this 
observation allows us to relate our finding for rigid surfaces to experi- 
mental results for water hydrating membranes 63 , reporting new types 
of water dynamics in thin interfacial layers, and water nanoconfined 
in different types of reverse micelles 64 , showing that the water 
dynamics is governed by the presence of the interface rather than 
the details (e.g., the presence charged groups) of the interface. 

In conclusion, we analyze the low-T phase diagram of a water 
monolayer confined between hydrophobic parallel walls of size L 
separated by h ~ 0.5 nm. We study water fluctuations associated 
to the thermodynamic response functions and their relations to the 
loci of TMD, TminD. For each response function we find two loci of 
extrema, one stronger at lower- T and one weaker and broader at 
higher- T. These loci converge toward a critical region where the 
fluctuations diverge in the thermodynamic limit, defining the 
LLCP. We calculate the Widom line departing from the LLCP based 
on its definition as the locus of maxima of C and show that it coin- 
cides with the locus of strong maxima of the response functions. We 
find that the LLCP belongs to the 2D Ising universality class for L — » 
°°, with strong finite-size effects for small L. Surprisingly, the finite- 
size effects induce the LLCP universality class to converge toward the 
bulk case (3D Ising universality class) already for a system with a very 
pronounced plane asymmetry, i.e. a water monolayer of height h ~ 
0.5 nm and L/h~ 50. For normal liquid, instead, this is expected only 
for much smaller relative values of L (L/h £ 5). We rationalize this 



result as a consequence of two properties of the HB network: (i) its 
high cooperativity, that enhances the fluctuations, and (ii) its low 
coordination number, that makes the first coordination shell for the 
monolayer and the bulk similar. 

Methods 

The model. We consider a monolayer formed by JV water molecules confined in a 
volume V = hi} between two hydrophobic flat surfaces separated by a distance h, with 

V/2V>Vo — 42 A , where v 0 is the water excluded volume. Each water molecule has 
four next-neighbours 7 . We partition the volume into JV equivalent cells of height 
/i~0.5nm and square section with size r= *J L 2 /N, equal to the average distance 
between water molecules. By coarse-graining the molecules distance from the 
surfaces, we reduce our monolayer representation to a 2D system. We use periodic 
boundary conditions parallel to the walls to reduce finite-size effects. We simulate 
constant JV, P, T, allowing V( T, P) to change, with each cell i = 1 , . . . , JV having number 
density Pj = p(T,P) = JV/V</? 0 = l/v 0 corresponding to a mass density — 1 g/ cm 3 . 
To each cell we associate a variable «, = 0 («,- = 1) depending if the cell i has pj/p 0 ^ 
0.5 (pj/po > 0.5). Hence, n, is a discretized density field replacing the water 
translational degrees of freedom. The water-water interaction is given by 

Jlf = Ufa) -JN HB -J a N coor . (9) 

9 

The first term, summed over all the water molecules /' and j at O-O distance r,p has 
U{r) = =o for r < r 0 = \J Va/h — 2.9 A (water van der Waals diameter), 
U(r) = Ac [(r 0 /r) 12 - (r 0 /r) 6 ] for r =: r a with c = 5.8 kj/mol, and U(r) = 0 for r> r c = 
25r 0 (cutoff). 

The second term represents the directional (covalent) component of the hydrogen 
bond (HB), with //4c = 0.5, N HB = njnjS^^ . number of HBs, with the sum 

over n.n., where cth= 1, ...,qis the bonding index of molecule / to the n.n. molecule;', 
with S a b = 1 if a = b, 0 otherwise. Each water molecule can form up to four HBs. We 
adopt a geometrical definition of the HB, based on the OOH angle and the OH — O 
distance. A HB breaks if OOH > 30°. Hence, only 1/6 of the entire range of values [0, 
360 '] for the OOH angle is associated to a bonded state. Therefore, we choose q = 6 to 
account correctly for the entropy variation due to the HB formation and breaking. 
Moreover, a HB breaks when the OH — O distance > r max — r Q H — 3. 14 A, where r Q H 
— 0.96 A and r max — 4.1 A. The value of r max is a consequence of our choice tij = 0 for 
pJp Q < 0.5, i.e. ij/2>r^, implying that = 0 when r« > r 0 v / 2 = 4.10 A = r max . 
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The third term of the Eq.(9) accounts for the HB cooperativity due to the quantum 
many-body interaction 65 , with / ff /4e = 0.05 andN COO p = y^ , n f ^ff*,^' where 

(/, k)i indicates each of the six different pairs of the four indices cTh of a molecule i. The 
value / ff <C/ is chosen in such a way to guarantee an asymmetry between the two 
components of the HB interaction. To the cooperative term is due the O-O-O 
correlation that locally leads the molecules toward an ordered configuration. In bulk 
water this term would lead to a tetrahedral structure at low P up to the second shell, as 
observed in the experiments 66 . An increase of Tor P partially disrupts the HB network 
and induces a more compact local structure, with smaller average volume per 
molecule. Therefore, for each HB we include an enthalpy increase Pv HB , where v HB /v 0 
= 0.5 is the average volume increase between high-p ices VI and VIII and low-/? 
(tetrahedral) ice Ih. Hence, the total volume is V= V 0 + Nhb^hb> where Vq ^ Nv 0 is a 
stochastic continuous variable changing with Monte Carlo (MC) acceptance rule 46 . 
Because the HBs do not affect the n.n. distance 66 , we ignore their negligible effect on 
the U(r) term. Finally, we model the water-wall interaction by excluded volume. 

The observables. The LLCP is identified by the mixed-field order parameter M and 
not by the magnetization of the Potts variables Gu as in normal Potts model. M is 
related to the configuration of the system by the relation 



M = 



+ 5 U(r) 



m ' Mi 



(10) 



where v = V 0 /N and s is the mixed-field parameter. M is therefore a linear 
combination of density and energy. 

Thermodynamic response functions are calculated from 



Kt — 



1 [S(V)\ _ (AV 2 ) 
(V)\dPj T k B TV 



(11) 













\v Gr)J 



— apKr - 



1 d 2 V 
VdPcT ' 



1 (BV\ fdV 
V^\dPj T {df 



dK T 



1 tfv 

VdTdP ' 



(16) 



Following 39,67 the line of extrema in density (TMD and TminD lines) is characterized 
by dp = 0, hence, dcc P = 0 along the TMD line. Therefore, 



0 = da P = 



dT-i 



<«-=| OT ] dT+(-— ) dP(17) 



where the index "ED" denotes that the derivatives are taken along the locus of 
extrema in density. So, the slope dP/dT of TMD is given by 



Tf 1 



BPdTI ^ 



(18) 



from which, using Eq. (15) with tx P = 0, we get Eq. (1). The Eq. (18) holds as long as 
both (ddp/dP) T and {8a P l8T) P do not vanish contemporary, as it occurs along the 
Widom line, where the loci of strong minima of Op overlap. For this reason the 
intersection between the Widom line and TminD line does not imply any change in 
the slope (dP/ 8T) TmiaD . 

To calculate Eq. (2) we start from Cp and otp written in terms of Gibbs free energy 



Q> 
T 



' BT 2 



V-j.,, 



(19) 



from which results 



and 



dT 



(AH 2 ) 



(12) 



as long as the volume and energy distributions are not clearly bimodal, i.e. excluding 
the values of T and P where the phase coexistence is observed, based on the definition 
of M. Here AO = O - (0), for O — V, H and, H is the enthalpy of the system. 

The Monte Carlo method. The system is equilibrated via Monte Carlo simulation 
with Wolff algorithm 46 , following an annealing procedure: starting with random 
initial condition at high T, the temperature is slowly decreased and the system is re- 
equilibrated and sampled with 10 4 + 10 5 independent configurations for each state 
point. The thermodynamic equilibrium is checked probing that the fluctuation- 
dissipation relations, Eq. (11) and (12), hold within the error bar. 



rP 



1 >, T\dPj T ~ 



OLp-V 



5 , 

— ( Vzp) 



W) 1 



h 



(20) 



(21) 



Substituting in Eq. (18) we get the Eq. (2) at the TMD. Moreover, because of ol p = 0 
at the TMD line, from the last equivalence of Eq. (20) we get 



Jf2 



(22) 



The histogram reweighting method. The probability Q N (M) is calculated in a 
continuous range of T and P across the c Max line. We consider an initial set of m e 
[10 : 20] independent simulations within a temperature range ATk B /(4c) ~ 10 ~ 4 and 
a pressure range APv 0 /(4c) ~ 10~ 3 . For each simulation i = 1, tn we calculate the 
histograms h,{u, p) in the energy density-density plane. The histograms h- t {u, p) 
provide an estimate of the equilibrium probability distribution for u and p\ this 
estimate becomes correct in the thermodynamic limit. For the NPT ensemble, the new 
histogram h(u, p, P', /?') for new values of /?' — l/k B T' and P' close the simulated 
ones, is given by the relation 57 



J2H i N je -l} i (u + P i /p)N-Q 



where iV, is the number of independent configurations of the run i. The constants C„ 
related to the Gibbs free energy value at T, and P,, are self-consistently calculated from 
the equation 57 



Q=-G(P ij ft)/fc B T. (14) 



We choose as initial set of parameters Q = 0. The parameters C, are recursively 
calculated by means of Eq. (13) and (14) until the difference between the values at 
iteration k and k + 1 is less then the desired numerical resolution (10~ 3 in our 
calculations). Once the new histogram is calculated, Q^(M) at T, and P f is calculated 
integrating h{u, p, P t , /?,) along a direction perpendicular to the line p + su. 

Thermodynamic relations. We report here the calculations for the thermodynamic 
relations in Eq. (1), (2), (3) and (4) 39 . To verify the relation (4) we calculate the 
derivative of K T along isobars 



— J_ (dV\ (dV\ 
~ V 2 KdTJpKfiPJT 



( <3t)p~ dT\ V {bp)t 

and the derivative of <x P along isotherms 



1 8 2 V 
V dPrT ~ 



-v-pKr-^jp^ (15) 



from which, using Eq. (18), we get the Eq. (3). 
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